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sipative lattices but the scaling exponents are different and seemingly nonuniversal. In contrast, 
backward LVs (obtained via Gram-Schmidt orthonormalizations) present approximately the same 
scaling exponent in all cases, suggesting it is an artificial exponent produced by the imposed or- 
thogonality of these vectors. We are able to compute characteristic LVs in large systems thanks to 
a 'bit reversible' algorithm, which completely obviates computer memory limitations. 
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I. INTRODUCTION 

In dynamical systems the sensitive dependence on ini- 
tial conditions is readily quantified by the Lyapunov ex- 
ponents (LEs) , which measure the average growth rate of 
infinitesimal perturbations 1]. Spatially extended sys- 
tems may exhibit spatiotemporal chaos (STC) and the 
spectrum of LEs is a good indicator of the extensivity 
with the system size [2, The directions in tangent 
space associated with the LEs are generically called Lya- 
punov vectors (LVs). These vectors convey important 
dynamical information. For example, LVs have been use- 
ful to discover and quantify the so-called hydrodynamic 
modes , to study extensivity properties @, Q and 
to address predictability questions in weather forecast- 
ing [j| [l(|, among other applications. 

For many dissipative one-dimensional models [ill, E3 
it is known that, after a suitable logarithmic transfor- 
mation, the infinitesimal perturbation associated with 
the largest LE, the main LV, belongs to the universality 
class of the stochastic Kardar-Parisi-Zhang (KPZ) equa- 
tion [l3j of surface growth. The mentioned logarithmic 
transformation associates a "surface" with the LV, lead- 
ing to many interesting consequences; for instance: the 
scaling of finite-size corrections and self averaging proper- 
ties of the LEs. The "surface picture" has demonstrated 
to be very powerful and has been used to analyze finite 
perturbations [I3 - [T6| and singular vectors [lTj in STC. 

The only homogeneous extended systems where the full 
correspondence between the main LV and the KPZ scal- 
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ing is known to break down are anharmonic Hamiltonian 
lattices. In Ref. [l8| it was determined, by numerical 
simulation of two different oscillator lattice models, that 
the main reason for the lack of KPZ scaling in Hamil- 
tonian systems can be traced back to the ubiquitous ex- 
istence of long-range spatiotemporal correlations in the 
observables that control the LV dynamics. In Ref. [l8[ 
the authors invoke the KPZ equation with a long-range- 
correlated noise (instead of white noise) as a minimal 
model that accounts for their observations. Neverthe- 
less, it remains unclear whether the KPZ equation with 
spatio-temporal long-range correlated noise is indeed the 
correct minimal model for the dynamics of the leading 
LV in one-dimensional Hamiltonian lattices. Further the- 
oretical progress is needed to clarify this issue. 

Recent studies [H, H(| have extended the analysis to 
LVs corresponding to the most unstable directions (not 
only the leading one) in several dissipative systems. As 
reasoned above, it is clear that Hamiltonian lattices de- 
serve a separate study due to the peculiar behavior al- 
ready observed for the main LV. 

We employ the so called [2l[ characteristic Lyapunov 
vectors (CLVs) proposed many years ago by Ruelle [22| 
because they reflect the bona-fide directions in tangent 
space (see below). CLVs have been recently employed 
to characterize several aspects of STC, such as spatio- 
temporal correlations and extensivity flfjl |20L [2311 , hyper- 
bolicity [U, [25| , and Oseledec splitting [26M28] . In ad- 
dition CLVs have nice properties that may support their 
use also for ensemble forecasting in atmospheric mod- 
els 03. 

The aim of this paper is to explore universality prop- 
erties (if any) of CLVs for Hamiltonian lattices, as these 
systems are already peculiar in what concerns the main 
LV. In addition, we present an algorithm, specially de- 
signed for Hamiltonian systems, to compute CLVs in 
large systems with modest computer resources. 
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This paper is organized as follows: Sec. |TT] describes the 
employed models and the relevant details of their numer- 
ical implementation. Sec. IIIII gives the relevant details of 
the "roughening surface" picture, as well as some results 
concerning the temporal evolution of the defined surface. 
In Sec. IIVI we investigate the spatial correlations of the 
CLVs. The discussion of the obtained results is made in 
Sec. El 



II. THE MODELS AND SIMULATION DETAILS 
A. Phase-space dynamics 

The reference Hamiltonian for the one-dimensional 
coupled anharmonic lattice models we are considering 
can be written as 



iV 



H = J2 



»=i l 



U(qi) 



(1) 



where N is the system size, and V(x) and U(x) are the 
nearest-neighbor interaction and on-site potentials, re- 
spectively. The particles are assumed to be of unit mass 
mi = 1. The phase space coordinates (displacement and 
momentum) are {qi,Pi}iLi'i periodic boundary conditions 
are assumed (q N+1 = q x )- In the following we shall con- 
sider two models: (i) the Fermi-Pasta-Ulam (FPU) (3 
model (29j ). characterized by V(x) = x 2 /2 + x 4 /4 and 
U[x) = 0, and (ii) the $ 4 model [30L l3lj|. characterized 
by harmonic interactions V(x) — x 2 /2 and by a double- 
well on-site potential U{x) — —x 2 /2 + x 4 /A. 

In our numerical simulations we have chosen as initial 
conditions the equilibrium value of the oscillators dis- 
placements, i.e. q(0) = {<Zi(0) = 0}£l 1: and momenta 
p(0) = {Pi(0)}iLi drawn from a Maxwell-Boltzmann dis- 
tribution at a temperature consistent with a given value 
of the energy density e = E/N. A value of e = 10 has 
been chosen for the FPU model, since it is known that 
its dynamics is strongly chaotic for e » 1 [30|, HH . For 
the $ 4 model we chose e = 5, as in Ref. (l8j . 



This linear evolution of infinitesimal perturbations im- 
plies the existence of a linear operator (resolvent or prop- 
agator) Ai that links perturbations at different times: 
ST{t) = M(t,t ) ■ 6T(t ). 

According to Oseledec's multiplicative ergodic the- 
orem [32j ] the remote past limit symmetric operator 
#b(t) = lim t0 ^_ oo [A4(t ! f ) -M*{t, t )] 1/[2( *- to)] exists 
for almost any initial condition T(to). The set of LEs 
is defined as X a = lnA Q , where {A Q } are the eigen- 
values of 3>b(£)- We label the LEs in decreasin g or der: 
Ai > A 2 > • • • > A 2 at. The standard procedure |33l. HiJ 
to compute the Af largest LEs resorts to periodic Gram- 
Schmidt-orthonormalizations of a set of Af offset vectors 
evolved by Eqs. ([2]). The time-averaged values of the log- 
arithms of the normalization factors yield the LEs {A Q }. 
The set of vectors right after each reorthonormalization 
{b Q (i)} are the eigenvectors of 3>b(i) [H[ and they are 
called backward LVs (BLVs) , following the nomenclature 
by Legras and Vautard [2lj]. Note that BLVs, apart from 
the main one bi(t), are not univocally defined because 
they depend on the scalar product for the orthogonal- 
ization (which also determines the adjoint operator A4*; 
Ai* = A4 T in Euclidean space). 

BLVs have the advantage of a straightforward calcu- 
lation as they are simply the byproduct of the standard 
method to compute LEs. However from the point of view 
of the physical meaning there is another set of vectors, 
the already mentioned CLVs (also known as covariant 
LVs), which univocally determine the direction in tan- 
gent space corresponding to each LE. The CLVs {g a (t)} 
were already defined by Ruelle in 1979 [H, Hl|. These 
vectors are independent of the definition of the scalar 
product and readily signal the intrinsic stable and unsta- 
ble directions. As a result CLVs are covariant with the 
linear dynamics, g Q (t) oc Ai(t, io)ga(io)> wherewith it is 
automatically guaranteed that the LEs are recovered in 
both, past and future, time limits: 



lim (t-t )- 1 \n\\M(t,t )g a (t )\\= \ a 



C. Important numerical issues 
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B. Tangent-space dynamics 

To study the local dynamical stability of our sys- 
tem we introduce the infinitesimal perturbations of the 
trajectory T (t) along all possible directions (position 
and momentum axes) of the phase space as ST(t) = 
(<5q (t) , 5p (t)), thus defining the 2iV-dimensional tangent 
space. These infinitesimal perturbations are governed by 
the linear equations 

d 2 H . d 2 H . 
oqi = » » oq, + -K-r°pi (2a) 
dqi dpi dpf 

■ d 2 H 8 2 H 

dqf dq l dp l 



The evolution of the phase space trajectory T(t) = 
(q(t),p(t)) is obtained integrating the 2N first-order 
Hamilton equations of motion. We have used a symmet- 
rical version of the velocity Verlet integrator specially 
suited for long-time simulations (3tJ , see Eq. jATTJ in the 
Appendix. The adopted time step value At = 0.01 as- 
sures a faithful representation of the Hamiltonian flow 
and a driftless average value of the total energy E with a 
fluctuation level of AE/E w 10 -3 — 10~ 4 depending on 
the system size. 

The computation of CLVs is not straightforward. We 
have used the method proposed by Wolfe and Samelson 
in Ref. [IH, wherein all relevant details are given. To 
find the A/th CLV one needs to compute: (i) the first Af 
BLVs and (ii) a set of Af — 1 vectors (forward LVs) in- 



tegrating backwards in time the perturbations that obey 
the adjoint operator of the linear dynamics (proceeding 
as in the case of BLVs, but using the transpose of the 
Jacobian matrix instead). 

The problem of the time-reversed integration is that, 
although the employed algorithm to integrate the equa- 
tions of motion derived from the Hamiltonian is explic- 
itly time reversible, the computed trajectories in phase 
space, obtained after reversing all momenta, do not coin- 
cide with those traced by the time-forward motion (due 
to the effect of round-off errors and chaos sensitivity). 
We have solved this problem by a suitable use of integer 
arithmetic operations, which suppresses unwanted nu- 
merical effects. Thus the original phase-space trajectory 
can be exactly traced back (see the Appendix for techni- 
cal details). Note that our procedure consumes almost no 
computer memory because it only requires to store the 
set {b Q (i)}^ 1 at the times where the A/th CLV is going 
to be computed. If the bit reversible algorithm were not 
used the state of the system would have to be recorded 
periodically to allow a faithful trajectory backtracking. 
Furthermore, if instead of the Benettin method the QR 
method were used (as in Ref. [131), the periodical storage 
of the R matrices required by the latter would quickly 
lead to a computer memory overflow. 

III. SURFACE GROWTH PICTURE 

For every CLV (likewise for BLVs) g a (t) — 
(<5q( Q ' (t) , (5p( Q ) (i)) it is convenient to define an asso- 
ciated "surface" 

h<t ] (t) = In y/[Sq^ (tW + lSp^ {t)Y (4) 

where, as before, index i = 1, ...,JV plays the role of 
space. Hereafter we refer to {h\ } as surfaces because a 
relation, for a = 1 , between this kind of log-transformed 
LV and stochastic surface growth equations was discov- 
ered in Refs. [ill EH (and proposed in Ref. [Hj] for Hamil- 
tonian lattices). Since U a \t) = (l/N)J2? =1 4 a \t) is 
the logarithm of a norm, the ath LE corresponds to 
the average velocity of the corresponding ath surface, 
(dh^ a \t)/dt) — A Q . Fig. Q] presents a specific example 
of the nice properties of CLVs. Perturbations at t = 
along the first five CLVs are let to evolve freely, i.e. obey- 
ing Eqs. ([2]), and the mean height of the associated sur- 
faces are computed versus time. The average velocities 
are fairly close to the corresponding LE values. 

The main LV gi = bi in spatio-temporal chaotic sys- 
tems is strongly localized in space [39l-|41j|. and transfor- 
mation (J4j) allows to unfold the spatial structure of the 
vector, which would be otherwise hidden close to zero. 
The localization of the main LV is dynamic, i.e. there is 
a slow wandering of the localization region. In Ref. [l8[ 
it was demonstrated by means of numerical simulations 
that, contrary to dissipative systems and other systems 
with STC, the surface associated with gi does not fall 
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FIG. 1: (Color online) Evolution of the mean height of sur- 
faces (h a (t)) associated with small perturbations initialized 
along characteristic LVs with indices a — 1, . . . , 5 from top to 
bottom (FPU model, N — 32 and energy density e = 10). The 
curves are averages over 100 realizations. The fitting slopes 
m a are in good agreement with the corresponding values of 



into the universality class of the KPZ equation. This fact 
is attributed to the presence of long-range correlations in 
space and time in Hamiltonian lattices (18(. 

For the FPU model we depict snapshots of surfaces 
corresponding to the first and second LVs in Figs. [2] (a) 
for the BLVs and (b) for CLVs. The vectors are strongly 

localized; notice that h\°^ , Eq. Q , is a logarithmic vari- 
able. In Figs. [Hc,d), we plot the time evolution of the 
localization sites corresponding to BLVs and CLVs, for 
a = 1, . . . , 10. We define the localization site as the posi- 
tion i where takes its largest value at a given time. It 
can be readily seen that the maxima corresponding to the 
BLV-surfaces are scattered all over the spatial domain, 
just as in the case of the coupled-map-lattice (CML) 
studied in Ref. [l9| where it was argued that this be- 
havior is a byproduct of the orthogonalization procedure 
and not a physical property of the perturbation dynam- 
ics. On the contrary, CLVs present much more correlated 
localization sites. 

The scaling properties of LVs corresponding to LEs 
smaller than the first one have been recently reported 
for spatio-temporally chaotic dissipative systems [ltl ■ 
These works revealed that LV-surfaces are piecewise 
copies of the main one. This is readily seen defining the 
difference-field 

fi* S hP-hP. (5) 

We show in Fig. [3l using CLVs, that this qualitative fea- 
ture holds for Hamiltonian lattices as well (it is irrelevant 
that the first LV does not belong to the KPZ universal- 
ity class). As in dissipative systems [l9l [20|. the typical 
plateau size of rf°^ decreases as a grows, and beyond 
some a m ax this simple picture does not hold: Figure [3fb) 
shows that for a = 8 the plateaus are smaller than for 
a = 2, Fig. Eta). 
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FIG. 2: (Color online). (a,b) Snapshot of the LV-surfaces 
corresponding to a = 1, 2 for the FPU lattice with N = 512. 
(c,d) Space-time evolution of the localization sites for the first 
10 LVs. The position of the localization site for the main LV 
is indicated by filled circles. Other symbols correspond to 
a = 2,..., 10.' 
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FIG. 3: (Color online) Snapshots of the field T]l at four 
different times for (a) backward LVs and (b) characteristic 
LVs. (c,d) Same as (a,b) but for the field T]!j. . Note that the 
typical length of the plateaus is smaller in this case. (FPU 
model with N = 512 and e = 10.) 



IV. SPATIAL STRUCTURE 



In this section we perform a quantitative descrip- 
tion of the spatial correlations of the LV-surfaces 
h[ a \t). We compute the stationary structure fac- 
tor S a {k) = lim t ^ 00 (^ Q )(/c,t)^( Q )(-fc,i)), where 
U a \k,t) = N- 1 / 2 J2 j exp(27nkj)h ( J a) (t), with (...) in- 
dicating an average over different system trajectories 



(which correspond to different random initial conditions) . 

Figures (Ha,b) and 03a, b) show the structure factors 
of a representative set of LV-surfaces for FPU and $ 4 
models, respectively. For the main LV surface, a = 1, 
the short wavenumber scaling exponent 7 of the struc- 
ture factor (S(k) ~ fc 7 ) is clearly different from the one 
expected for the KPZ universality class in one dimension 
(7 = —2) and the obtained values, 7 = —2.5 (FPU) and 
7 = —2.6 ($ 4 ), are consistent with those reported by 
Pikovsky and Politi |18| . 

As explained in Sec. Ill Bl for a > 1 one must dis- 
tinguish between backward and characteristic LVs. Fig- 
ures 0Ja,b) an d [3a, b) evidence that both vector types 
indeed have different spatial structures. The structure 
factors of BLV-surfaces asymptotically decay with expo- 
nents 7 = -0.9 (FPU) and 7 = -1.1 ($ 4 ), which are 
close to the value 7 = —1 reported for dissipative sys- 
tems in Refs. [l|| These results suggest that the 
value 7 — 1 for BLVs in our Hamiltonian systems has 
a geometric origin and is related to the Gram-Schmidt 
orthonormalization. However, CLVs display exponents 
7 = -1.4 (FPU) and 7 = -1.3 ($ 4 ) which, at least for 
the FPU model, are different from the values —1.2 or 
— 1.15 reported in dissipative systems [i~9l [20j. 

The fact that LVs corresponding to the most expand- 
ing directions are (in the surface representation) piece- 
wise copies of the main LV translates into the existence 
of crossover wavenumbers where the structure factors 
bend. Each structure factor S a presents a knee at a 
certain wavenumber fc* that is related to the typical 
plateau length of the difference field rf a \ For both, 
FPU and $ 4 , models fc* scales with a approximately 
as fc* ~ [{a - 1/2) /N]°, with 6 » 1 (as in Ref. [H for 
a dissipative CML). We choose to use a — 1/2 instead of 
a, as it is customary when plotting the Lyapunov spec- 
trum (this is actually irrelevant in the thermodynamic 
limit N — >• 00). This is verified in Figs.@|c,d) and[5Jc,d) 
through a data collapse of the structure factors via the 
scaling relation 

S a (fc)fc- 7 =.g(fc/fc*), (6) 
where g(u) =const for m>1 and g(u) ~ u a for uCl. 



V. DISCUSSION AND CONCLUSIONS 

In a recent paper [2(| we established (as suggested 
by [Hj]) a minimal stochastic model for LVs w(x,t) in 
spatiotemporal chaos: 



d t w(x, t) = £(x, t)w(x, t) + d xx w(x, t), 



(7) 



where £(x,t) is a white noise term that accounts for 
the chaotic fluctuations. The generic asymptotic solu- 
tion of Eq. ([7]) models the main LV, whereas "saddle 
solutions" (see 20] for details) correspond to CLVs with 
2 < a < ct max . This correspondence between a (linear) 
stochastic equation and LVs is supported by the existence 
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FIG. 4: (Color online) Structure factors S a (k) of the FPU 
model with TV = 4096, for (a) backward and (b) characteristic 
LV-surfaces, a — 1,4,8,16,32,64,120 from top to bottom. 
(c,d) Data collapse through scaling relation with 7 = —2.5 
and fitting parameters 9 = 0.9 (c) and 9 — 0.8 (d); index a 
increases from right to left. 
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FIG. 5: (Color online) Structure factors S a (k) of the $ 4 model 
with TV = 2048, for (a) backward and (b) characteristic LV- 
surfaces, a = 1,5,10,20,35,60 from top to bottom. (c,d) 
Data collapse through scaling relation ((6| with 7 = —2.6 and 
fitting parameters 9 — 1.2 (c) and 9 = 1.0 (d); index a in- 
creases from right to left. 



these multipliers are functions of the displacements {qi} 
with specific properties for FPU and $ 4 models. The 
minimal model proposed in [l8[ for the main LV-surface 
was the KPZ equation [or ([7]) reversing the Hopf-Cole 
transformation] with long-range correlated noise. Unfor- 
tunately, theoretical expressions for 7 [52, EH consider 
the KPZ equation with either spatially or temporally- 
correlated noise (but not both). The exponent 7 may 
take the same value with different combinations of spa- 
tial and temporal long-range correlations. Hence there is 
not a univocal relation between correlations and 7. 

Minimal models are important as they are more 
amenable to theoretical analysis, which should allow to 
distinguish different universality classes 1 . Our present 
work has the value of giving more constraints to the min- 
imal stochastic model for perturbations in Hamiltonian 
lattices. A minimal model should reproduce the scaling 
properties of both the main LV and sub-dominant LVs. 
These sub-dominant (a > 2) LVs have been the subject 
of the present study. A possible alternative to J7J pointed 
out in [i~2| is the time-reversible equation: 



d tt w(x, t) = £{x, t)w{x, t) + d xx w(x, t). 



(8) 



Future work is need to find out the true minimal model 
for LVs in Hamiltonian lattices. In any case, our results 
should be a guidance in the search of such a minimal 
model for Hamiltonian systems, since any suitable mini- 
mal model must produce surfaces with the scaling prop- 
erties in Figs. 0] and [5] 

Finally, a bit reversible algorithm, which operates with 
integer arithmetic, has been implemented to take ad- 
vantage of the time reversibility of Hamiltonian systems. 
Although trajectory reversibility is not strictly required 
to compute the CLVs (in previous works [l!| H(| CLVs 
have been computed in non-reversible dissipative sys- 
tems), our methodology makes use of the aforementioned 
reversibility to completely bypass the need to store the 
phase-space trajectory in order to maximize the efficiency 
of the computation of the CLVs from the intersection of 
backward and forward LV subspaces if only limited com- 
puter resources are available. 
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Appendix: Bit reversible algorithm 

Since the so called 'bit reversible' technique has been 
so far implemented using the standard Verlet integra- 
tor (which does not considers the momenta explicitly) 
and employed mainly for studies of time reversibility in 
Lennard- Jones fluids [45H47I ]. (although it has also been 
applied in cases in which the so-called smoothed-particle 
continuum mechanics becomes isomorphic to molecular 
dynamics, see Rcf. 48.]) we will give a brief explanation 
of its current implementation in order to make this pa- 
per self-contained. Starting from the initial condition 
r(0), the phase-space point r(Ai) is obtained by means 
of the symmetrical velocity Verlet integrator [33] written, 
in floating-point arithmetic, as 

p*(t) = Pi(°) + ^[*(o)]. ( A - la ) 

9i (At) = ft (0) + ^-Pi (^f) , (A.lb) 

Pi (At) = Pi (^pj + ^Fi[ qi (At)], (A.lc) 

where F{ is the total force on the ith oscillator. The 
bit reversible version of algorithm (|A.1[) employs an in- 
teger representation of phase space instead of the con- 
ventional continuous phase space. To accomplish such 
transformation, for the considered lattices the minimum 
distance by which the phase space is discretized is defined 
as AL = N/2 n , where N is the system size and 2 n is the 
largest integer value if n-bit integers are employed. Be- 
cause of the discretization, the phase space coordinates 
are represented by integers, i.e. {iq, ip}. Therefore the 



At 



iQi(At) 
iPi(At) 



1+1 ( At 1 
ip,(0)+ £ T ^fe(0)] (A.2a) 

7=2—1 *• ' Integer 



.At (At 



(A.2b) 



Into 



At 



'+ 1 r \ f -] 

]T — ^fe(At)] (A.2c) 

j—i— 1 ^ ^ Integer 



where F^ is a partial force from the j-th nearest neighbor 
on the i-th oscillator. It should be noted that the discrete 
coordinates {iq(t), ip(f)} are integers, and the actual 
phase-space coordinates are obtained as q(t) = iq(t)A£ 
and p(t) = ip(t)AL. The second terms in the right hand 
side of Eqs. (|A.2j) are calculated based on the continuous 
phase space variables {q, p} and the values in brackets 
{■jinteger are converted to integers; in this way the total 
momentum is exactly zero at all times during the simu- 
lation |45|. 

With the aforementioned implementation time re- 
versibility is achieved exactly throughout the simulations 
performed, which were rather long. To give an example, 
for the $ 4 model a simulation of 5 x 10 8 time steps, after 
a transient of 1.5 x 10 8 , was needed to obtain the re- 
ported results. The situation is definitely better for the 
FPU model, where only 5 x 10 6 time steps, with a tran- 
sient of 5 x 10 5 , were sufficient for the employed system 
size. Due to the conversion process to integer, energy 
is not exactly conserved. Nevertheless, the exact-time- 
reversibility of the integration algorithm precludes any 
systematic drift in the total energy. We have confirmed 
that the fluctuation of the total energy by the bit re- 
versible simulations is equivalent to that by conventional 
floating-point simulations. 
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